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ABSTRACT 

This paper examines the relative impacts on grid-averaged longwave flux transmittance (emittance) for marine 
boundary layer (MBL) cloud fields arising from horizontal variability of optical depth rand cloud sides. First, 
using fields of Land sat-inferred r and a Monte Carlo photon transport algorithm, it is demonstrated that mean 
all-sky transmitiances for 3D variable MBL clouds can be computed accurately by the conventional method of 
linearly weighting clear and cloudy transmittances by their respective sky fractions. Then, the approximations 
of decoupling cloud and radiative properties and assuming independent columns are shown to be adequate for 
computation of mean flux transmittance. 

Since real clouds have nonzero geometric thicknesses, cloud fractions A, presented to isotropic beams usually 
exceed the more familiar vertically projected cloud fractions A, . It is shown, however, that when A t ^ 0.9, biases 
for all-sky transmittance stemming from use of A, as opposed to A ( are roughly 2-5 times smaller than, and 
opposite in sign to, biases due to neglect of horizontal variability of r. By neglecting variable r, all-sky trans- 
mittances are underestimated often by more than 0.1 for A, near 0.75 and this translates into relative errors that 
can exceed 40% (corresponding errors for all-sky emittance are about 20% for most values of A, ). Thus, priority 
should be given to development of general circulation model (GCM) parameterizations that account for the 
effects of horizontal variations in unresolved r, effects of cloud sides are of secondary importance. 

On this note, an efficient stochastic model for computing grid-averaged cloudy-sky flux transmittances is 
furnished that assumes that distributions of r, for regions comparable in size to GCM grid cells, can be described 
adequately by gamma distribution functions. While the plane-parallel, homogeneous model underestimates cloud 
transmittance by about an order of magnitude when 3D variable cloud transmittances are £ 0.2 and by —20% 
to 100% otherwise, the stochastic model reduces these biases often by more than 80%. 


1. Introduction 

While there is a wealth of studies aimed at shortwave 
radiative transfer for inhomogenous clouds (Welch and 
Wielicki 1984; Davis et al. 1990; Barker 1992; Cahalan 
et al. 1994a), there is a dearth of studies regarding the 
impact of inhomogeneities on longwave (LW) radiative 
transfer (e.g., Harahvardhan et al. 1981; Ellingson 1982; 
Evans 1993). In remote sensing and climate modeling 
studies, clouds are usually assumed to be horizontally 
homogeneous and occasionally assumed to be black in 
the LW portion of the spectrum. For (60 km) 2 regions 
of marine boundary layer (MBL) clouds, however. Bar- 
ker et al. (1996) showed that small values of cloud op- 
tical depth r are often very abundant even when mean 
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t is much larger than 0 (see also Wielicki and Parker 
1992, 1994). Thus, since LW radiative transfer is an 
extremely nonlinear process for small r, it is reasonable 
to expect that grid-averaged LW transmittances and em- 
issivities, as required by general circulation models 
(GCMs), may be sensitive to unresolved horizontal in- 
homogeneities of cloud. 

A case in point supporting this expectation is Wielicki 
and Parker’s (1992) assertion that by neglecting trans- 
lucent clouds, the spatial coherence method (Coakley 
and Bretherton 1982) often underestimates MBL cloud 
fraction by 0.15-0.2. Heeding this, Luo et al. (1994) 
extended the spatial coherence method and deduced that 
for (250 km) 2 regions of marine stratocumulus clouds 
off the coast of South America, grid-averaged 11 -/xm 
emissivities for clouds only are often closer to the value 
of cloud fraction than to unity. For this to be true, not 
only must there often be substantial amounts of thin 
cloud amid readily observable thicker cloud, but the 
inhomogeneous nature of these clouds is likely impor- 
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Fig. t. Schematic 2D cross section of a cloud field (consisting of 15 cells) with geometric 
thicknesses h, and vertically projected cloud fraction 0.8. For lines-of-sight 30° from the zenith 
(fi ** 0.86), however, apparent cloud fraction is almost 100%. Note that distributions of optical 
depth encountered by zenith traveling radiances (e.g., those depicted on the left) are composed of 
Tj only, whereas those for off-zenith radiances (e.g., the one grazing cell 7) consist, in part, of ar, 
where a ^ 1. Also for off-zenith radiances, distributions of optical depth consist of spatially 
averaged values such as that associated with the radiance spanning cells 12 through 15. 


tant for LW radiative transfer and thus cloud heating 
rates. 

In addition to unresolved fluctuations in r, numerous 
studies have demonstrated the importance of cloud sides 
in radiative transfer for broken clouds fields (e.g., 
Harshvardhan et al. 1981; Ellingson 1982; Zuev and 
Titov 1995). Yet for LW radiative transfer, the relative 
impacts of cloud sides and horizontal fluctuations in r 
have not been delineated. 

The main purpose of this paper is to examine the 
impacts on grid-averaged LW transmittance (emittance) 
for inhomogeneous MBL clouds arising from cloud fluc- 
tuations at scales unresolved by GCMs. Section 2 gives 
a background discussion of LW flux transmittance for 
inhomogeneous MBL clouds. In section 3, optical 
depths inferred from Landsat data are presented briefly. 
These data are used in section 4 as input for a Monte 
Carlo photon transport algorithm, the results of which 
culminate in an assessment of the relative impact of 
cloud sides and horizontal variability of r. Section 5 
presents an approximate technique for computing grid- 
averaged cloudy-sky transmittance that assumes that r 
are distributed according to the gamma distribution 
function and that the independent pixel approximation 
(see Cahalan et al. 1994a,b; Barker 1996) applies for 
LW radiation. Concluding remarks are in section 6. 

2. Longwave transmittance for horizontally 

inhomogeneous cloud: A conceptual model 

To begin, assume that the intensity of a pencil of LW 
radiation is extinguished by clouds in accordance with 
the Beer-Bouger-Lambert law and that multiple scat- 
tering by droplets can be neglected since droplet single 
scattering albedo a> 0 is small (<0.5) and asymmetry 
parameter g is large (>0.9) for much of the atmospheric 
window. Thus, throughout this study, r symbolizes LW 
absorption optical depth, which equals (1 - <o 0 )T cxt 


where r exl is extinction optical depth. Furthermore, ex- 
tinction by gases is neglected in order to simplify the 
presentation and since only boundary layer clouds are 
considered; all clouds are assumed to be isothermal. 

Figure 1 shows a schematic diagram of the main con- 
cerns involving LW radiative transfer through horizon- 
tally inhomogeneous boundary layer clouds. Consider 
viewing this cloud (field) at different zenith angles, 0 
(/a = cos0). For simplicity, all quantities throughout this 
study are assumed to be azimuthal averages. The prob- 
ability of a line-of-sight being intercepted by cloud is 
a function of /a, minimized for jla = 1 and increasing 
monotonically as /l i decreases. This probability can be 
thought of as the zenith-angle-dependent cloud fraction 
A c {p). The form of Afp) depends on several factors 
including distributions of cloud size (Wielicki and 
Welch 1986), aspect ratio (Plank 1969), and spacing 
(Cahalan 1991). For nonovercast clouds, the strongest 
and weakest dependencies of A c (/a) on ju are probably 
associated with fields of towering clouds and MBL 
clouds, respectively. Given that observations and radi- 
ative fluxes are affected by A c (/a), when cloud fractions 
are reported for observations and GCMs, it is unclear 
what is being, and what should be, discussed. It seems 
likely that the most common interpretation of the term 
cloud fraction is the vertically projected value A c (l). 
This quantity, however, is likely neither that reported in 
cloud atlases nor that most meaningful for computation 
of radiative fluxes in GCMs (i.e., ID column models). 

Next, consider probability distributions of cloud op- 
tical depth r (normalized to the vertical as usual) for 
lines of sight along given zenith angles. Denote these 
distributions of r conditional upon p as p(r I /a). As can 
be inferred from the idealized clouds in Fig. 1, one can 
expect p(r \ 1) to be relatively broad but as /a — » 0, 
p(r\ jjl) will tend to become narrower and more sym- 
metric about the zenith-angle-dependent mean cloud op- 
tical depth r(/Lt). This is because averaging optical depth 
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along oblique lines through a cloud field is a form of 
horizontal smoothing that is somewhat analogous to the 
narrowing of p(r 11) as the resolution of a cloud field 
degrades (e.g., Schertzer and Lovejoy 1987). Likewise, 
the form of r(p) will depend on the geometry of the 
cloud field. For example, the clouds in Fig. 1 have thin 
wedges near their tops and sides that are exposed to 
lines of sight with p< 1. Thus, one would expect r(p), 
and the variance of r, to decrease slightly with decreas- 
ing p. 

Taking the cloudless and cloudy parts of a region 
together, the grid-averaged, all-sky distribution of op- 
tical depth conditional upon p is 

^P(r\p) = [1 - A c (p)]8(r) 4 A c {p)p(r\p), (1) 

where 5(r) is the integrand of the Dirac function. As- 
suming that 

p(r \p) dr = 1 and rp(r\p) dr - ~T{p), 

Jo Jo 

(2a) 


2>(t I /a) is also normalized as 

I !P(r \ p) dr = 1 — A c (p) 4 A c (p) = 1, (2b) 

Jo 

with grid-averaged zenith-angle-dependent optical 
depth of 


f 


t(P(t I p) dr 


= 04 Afp) 

Jo 


rp(r\p) dr = A c (p)T(p). (2c) 


Similar expressions exist, of course, for all moments of 
^(rl p). 

Letting radiance transmittances through cloud and 
clear air be e~ T/fX and 1, respectively, grid-averaged all- 
sky transmittance T for an isotropic distribution of in- 
cident radiation is 


-/7 


{[1 - A r (/0]S(r) 

+ A c (p)p(r I p)e~ Tl,i }p dp dr 


= l - A, + 2 \ I A c (p)p(rlp)e Tl>t p dp, dr, 
Jo Jo 


cloudy-sky contribution 


(3) 


in which 


-4 


A c (p)^ dp 


(4) 


could be referred to as the hemispherical cloud fraction. 


Like A£ 1), A c is almost certainly neither cloud fraction 
estimated by, and used in, GCMs nor cloud fraction 
reported in cloud atlases. Moreover, the cloudy-sky con- 
tribution in (3) shows that, in principle, radiation fields 
and cloud fraction cannot be decoupled in a straight- 
forward manner, even when scattering is ignored. Ste- 
phens (1988) discussed this issue and showed that for 
LW radiative transfer through opaque clouds, T can be 
expressed as a linear combination of clear and cloudy 
transmittances weighted by suitable clear and cloudy 
fractions. This approximation is ubiquitous to climate 
studies and, when applied to (3), it becomes 

T = 1 - A c 4 A t T cld , (5a) 

where 

T c id = 2 j j p{r\p)e r/ »p dp dr = 1 - e cld 

Jo Jo 

(5b) 

is grid-averaged cloud transmittance and £ cld is corre- 
sponding emissivity. 

For a plane-parallel, homogeneous (PPH) cloud of 
optical depth r, p(r\ p) = S(r - r), which, when sub- 
stituted into (5b), gives the familiar result 

T^ d = 2 j e~ T '»p dp = 2£ 3 (t) = 1 - £&, (6) 

Jo 

where E 3 (t) is the third-order exponential integral 
(Charlock and Herman 1976). Often, Tjf d is expressed 
as where <D is the diffusivity factor (Elsasser 1942; 
Stephens 1978). Quanhua and Schmetz (1987) ex- 
pressed D as a function of r such that e ~' pl is equivalent 
to (6). Wielicki and Parker (1994) and Barker et al. 
(1996) showed, however, that for MBL clouds, it is often 
the case that large fractions of area have small r, even 
for r » 0. Thus, at times, it can be expected that the 
nonlinear effects of averaging e T/fl or E 3 (r) over all r 
will yield results that differ greatly from T%(r). Given 
that cloud fractions reported in cloud climatologies in- 
clude thin clouds (Rossow and Schiffer 1991), and that 
future climatologies will report even thinner clouds (see 
Wylie et al. 1994), GCMs will have to account for these 
thin clouds in order to make validation of simultaneous- 
ly predicted cloud fraction and radiative budgets as un- 
ambiguous as possible. It stands to reason, therefore, 
that the radiative impact of all clouds should also be 
addressed by GCMs. Hence, it may be essential, at 
times, to utilize forms of r cld and A c that are less trivial 
than 7Pf d and A c { 1). 

Thus, the primary objectives of this study, in the con- 
text of MBL clouds only, may now be stated clearly as 
(i) to establish the applicability of (5a), (ii) to deduce 
the necessity for parameterizing A c , and (iii) to establish 
a suitable parameterization for r dd . Investigations are 
conducted using fields of Landsat- inferred cloud optical 
depths as input to a Monte Carlo photon transport al- 
gorithm. 
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Scene B2 Scene B1 3 



1 1 .5 pm Absorption Optical Depth 



0 0.25 0.5 1 2.5 5 10 25 50 


Fig. 2. Absorption optical depths r for 1 1.5-/Lim radiation for four scenes used in this study. 
These values were transformed from Landsat-inferred 0.83-/xm optical depths by multiplying by 
0.46. Table 1 lists some information about these images. 


3. Data 

Fields of optical depth inferred from 45 Landsat im- 
ages of MBL clouds were employed [see Barker et al. 
(1996) for a summary]. Each image is 60 km 2 , of which 
41 consist of 2048 2 pixels while the others consist of 
1024 2 pixels. They were presented originally by Harsh- 
vardhan et al. (1994), who used 0.83-/zm nadir radiances 
to derive cloud extinction optical depths t 083 at hori- 
zontal resolution of either 28.5 or 57 m (Wielicki and 
Parker 1994). Thus, each image has its own p{r \ 1). Use 
of these p(r 1 1) for radiative flux calculations seems ad- 
equate for at least two reasons. First, at these resolutions, 
the vast majority of individual clouds are resolved very 
well (Wielicki and Welch 1986). Second, since the am- 
plitude of variations in r are known to decay rapidly 
for spatial scales less than —500 m (e.g., Cahalan and 


Snider 1989), fluctuations at scales less than —60 m are 
likely to be inconsequential for radiative transfer cal- 
culations. 

Assuming the effective radius of cloud droplets r e to 
be 10 gm (Han et al. 1994), extinction optical depth for 
wavelength 1 1.5 pm is approximately equal to 0.78r O83 , 
and for 1 1.5-/im radiation is approximately 0.41 (Hu 
and Stamnes 1993). Hence, for this study, values of r 0 83 
are transformed into absorption optical depths for 1 1.5- 
pm radiation as 

r - (1 - 0.41)(0.78 )t O83 - 0.46t 083 . (7) 

While results are presented for all 45 scenes, addi- 
tional details are provided for the four scenes shown in 
Fig. 2: two examples each of broken stratocumulus and 
scattered cumulus. Table 1 lists information about these 





15 December 1997 


BARKER AND WIELICKI 


2789 


Table l. Summary of the four Landsat scenes shown in Fig. 2. 
Here A t is vertically projected cloud fraction, t is mean 11.5 
absorption optical depth for clouds only, and v is the gamma function 
parameter obtained from Landsat-inferred values of r: v = (f/cr) 2 , 
where it is standard deviation of r for clouds only. These quantities 
correspond to zenith radiances (i.e., /x = 1). 


Scene 

Date 

(d/mo/yr) 

Lat/Long 

A, 

f 

V 

B2 

10/7/87 

31.17°N/129.45°W 

0.644 

1.58 

1.251 

B 13 

13/7/87 

20.71°S/74.80°W 

0.974 

3.19 

1.068 

C7 

8/7/90 

33. 18°N/33.81°W 

0.256 

3.98 

0.669 

C12 

15/6/87 

7.22°S/1 15.58°W 

0.226 

2.79 

0.189 


scenes. Detailed results are not presented for an overcast 
example because their distributions of r are sufficiently 
narrow and r are sufficiently large (Barker et al. 1996) 
that LW transmittance and emissivity biases arising 
from the PPH assumptions are minor. 


4. Monte Carlo experiments 

This section presents results from a 3D Monte Carlo 
(MC) photon transport algorithm (Barker and Liu 1995) 
initialized with fields of Landsat-inferred r (the 11.5- 
p m absorption optical depths). In the MC experiments, 
cyclic horizontal boundary conditions were assumed 
and since absorption optical depths were used, a> 0 = 0. 
Cloudy pixels were modeled as vertically homogeneous 
columns with constant top elevation and variable geo- 
metric depth prescribed (in meters) as 

h = 75.9r 2/ \ (8) 


which is a good approximation to the Minnis et al. 
(1992) curve fit. While this prescription for cloud thick- 
ness is certainly not perfect, it is shown later that it 
yields reasonable standard deviations for h. All-sky MC 
transmittances T were computed by showering the ar- 
rays with isotropic distributions of ~3.7 million pho- 
tons. From (3), T can be defined as 


T = 


A,. + 



A c (p)T(p)p dp , 


(9) 


where T(/a) is zenith-angle-dependent, mean cloud 
transmittance. Since r was accumulated for each photon, 
A t . is defined as the fraction of photons that accumulated 
r > 0. To generate the ^.-dependent functions, /a- specific 
simulations were performed at select angles for /a be- 
tween 0.2 and 1.0. 

First, consider the approximation of decoupling cloud 
fraction and cloud transmittance. This leads to (5) and 
was tested here by simply assessing how well 


[ A c (fi)T(fi)fjL dfi = A, 

2 | dp. 

Jo 

■>« 

= A TW 
— A c eld * 



Figure 3 shows the left- and right-hand sides of (10) 



Cloudy-sky contribution in MC model 

Fig. 3. Cloudy-sky contribution to overall transmittance as com- 
puted by the Monte Carlo algorithm (abscissa) plotted against the 
product of Monte Carlo determined hemispherical cloud fraction A t 
and mean cloudy-sky transmittance Tffi (ordinate) for 45 Landsat 
scenes [see Eq.( 10)]. Points on the diagonal line indicate near equality 
in (10). 


plotted against each other for all 45 scenes. Clearly, the 
error in this decoupling approximation is very small and 
leads to only a slight, but systematic, overestimation of 
the cloud contribution to overall transmittance (<0.01 
at most). 1 That the largest differences tend to be asso- 
ciated with the smallest values of A r (I) (not shown) 
indicates further that, for all-sky transmittances, this ap- 
proximation is adequate. This, therefore, is taken as jus- 
tification to use (5). 

The next two tests have a bearing on how well p(r\ p) 
can be approximated by simply p(r II). The fitted re- 
lations presented below were confined to /a > 0.2: they 
tended to break down often for p < 0.2. This limitation 
poses little problem for flux quantities, however, as 96% 
of an isotropic beam is within p ^ 0.2. 

Figure 4a shows f{p) for the four scenes in Fig. 2. 
The value of f(p) tends to decrease slightly as p de- 
creases on account of increased exposure of thin corners 
on the rectangular columns of cloud. In most cases, r(p) 
can be fit very well with the regression line 

J(p) ~ r(l) + 0,(1 - /a), (11) 

where a x is a coefficient determined by least-squares 
linear regression (as are all coefficients in this section). 
Figure 4b shows that for the most part, a ] 555 0 can be 
expected. Moreover, scenes that comply worst with (11) 
(i.e., smallest coefficients of determination R 2 ) are as- 
sociated with very small, and irrelevant, values of a x . 
The outlier with A,(l) ~ 0.53 was observed at a solar 
zenith angle 0 o of 69° (most others were at 0 o ~ 30°) 


1 If A(/x) = up 1 , which is often an excellent fit (as shown later), 
and T(/a) = a + + y yx 2 , which is also a good representation, it 

is straightforward to show that 2 /,’> A t (pYT(fi) ijl dfjL < A r 7^ if b < 
0. which is true for the clouds considered here. 
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T(J1)=T(1)+S 1 (1 -H) 



4 , 0 ) 


Fig. 4. (a) Solid lines are r {/jl) as determined by the Monte Carlo 
algorithm for the four scenes shown in Fig. 2. Broken lines are least- 
square linear regression fits of the form listed at the top of (b), [see 
also (11)]. (b) Slope coefficients a, for the fits of (11) to 45 scenes 
as a function of corresponding A,(l). 


V(H)=V(1)H*2 



Fig. 5. (a) Solid lines are as determined by the Monte Carlo 
algorithm for the four scenes shown in Fig. 2. Broken lines are least- 
square linear regression fits of the form listed at the top of (b) [see 
also (12)]. (b) Exponent coefficients a 2 for the fits of (12) to 45 scenes 
as a function of corresponding A,(l). 


and this could be problematic (cf. Loeb and Coakley 
1997). Roughly speaking, the quantity of concern for 
flux transmittances is similar to /xe*^ and so minor 
changes in t (/x) with respect to fx are negligible because 
jxe^^^ generally approaches zero rapidly as /x decreas- 
es. 

Let the quantity v{fi) = [t(/x)/ct(/x)] 2 , where (^(fx) 
is zenith-angle-dependent variance of r, be a measure 
of relative magnitude of horizontal variability. Figure 
5a shows that v(/x) tends to increase slightly with de- 
creasing fx. Since r(j lx) ~ const, this means that (^(/x) 
decreases as /x decreases, which is not surprising and 
follows from the discussion in section 2 regarding 
smoothing via horizontal sampling by off-zenith radi- 


ances. The curves for v(/x) in Fig. 5a are described well 
by 

v(jx) ^ p(1)/x% (12) 

where a 2 is again a regression coefficient. Figure 5b 
shows a 2 as a function of A c ( 1) for the 45 scenes. In 
general, the smaller A ( .(l), the more smoothing takes 
place for off-zenith trajectories. Again, (12) fits worst 
when a 2 is small, implying that (12) is an adequate 
model that tends to break down only when the trend it 
attempts to capture is of negligible importance. As in 
Fig. 4b, the outlier with A c (l) ~ 0.53 is again an outlier 
in Fig. 5b. The value of a 2 for scenes B13 and Cl 2 are 
relatively large at about -0.3, and from Fig. 5a it can 
be seen that this indicates that for most scenes and most 
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Fig. 6. Solid lines are A {ix) as determined by the Monte Carlo 
algorithm for the four scenes shown in Fig. 2. Broken lines are least 
square linear regression fits of the form shown in (13). 


1 1 , v(fi) ~ i'(l): only at small fx does a{fx) become 
significantly small, but the contribution of the associated 
transmittances to the hemispherically integrated trans- 
mittance is quite small also. Hence, there appears to be 
very little reason to encumber a parameterization with 
v(jjL). Furthermore, the jx dependence of r(fx) and a^ifx) 
tend often to be in opposition: as fx decreases, reducing 
f(fx) enhances transmittance, while enhancing v{jx) re- 
duces it (cf. Barker et al. 1996). 

For simplicity, \x dependencies of r(p.) and v(fx) are 
neglected hereinafter and referred to as just t and v } 
which are taken to be equivalent to r(l) and v(\). Thus 
far, the results of this section indicate, fortunately, that 
simple parameterizations of 7^, such as by the con- 
ventional independent pixel approximation (IPA) (Ca- 
halan et al. 1994a, b; Barker 1996), can be applied with 
confidence. 

The next stage examines the necessity of having to 
use A t . as cpposed to, for example, simply A <; ( 1 ), which 
may be the cloud fraction most people think of as pre- 
dicted by GCMs and reported in cloud climatologies. 
Figure 6 shows values of A c (fx) for the scenes shown 
in Fig. 2. It also shows that for most A c (/x) are ap- 
proximated very well by 

A v {fi) ~ A t (Dti a \ (13) 

where a 3 is a regression coefficient that depends on dis- 
tributions of cloud size, spacing, and aspect ratio. Note 
that the magnitude of a 3 for scene C12 is the largest of 
the four shown, and from Fig. 2 it can be seen that its 
clouds are small and spread quite uniformly over the 
field. Thus, it is easy to see why for decreasing fx, A r .(/ x) 
increases so much relative to A c (l). Clouds in scenes 
B2 and C7, however, are more clustered and this results 
in smaller values of a y (clear lines-of-sight are difficult 
to close until very small fi). Figure 7a shows that for 
the 45 scenes, is correlated fairly well with vertically 


/4 c 0l)=/4 c (1)g*3 


CO 




Fig. 7. (a) Exponent coefficients a, for the fits of (13) (see Fig, 6) 
to 45 scenes as a function of corresponding A,(l). Broken line is a 
least-square linear regression fit defined in ( 15). (b) Coefficients of 
determination R 2 for the fits listed at the top of the plot as a function 
of A,( 1 ). 


projected cloud fraction A,(l). For reasons just alluded 
to, the tendency in Fig. 7a is the smaller A t (l), the 
greater the dependence of A r (A0 on /*• Figure 7b shows 
R 2 that result from fitting (13) for all 45 scenes as a 
function of A ( (l). Most cases exhibit R 1 > 0.98, im- 
plying excellent fits. As with v(fx) and r(fx), the few 
cases with poor fits are near overcast with very small, 
and irrelevant, values of a y Again, the scene with A,(l) 
~ 0.53 has an anomalously small value of a y . Oddly, 
however, Fig. 7b shows that for this scene the fit in (13) 
is excellent. 

Substitution of (13) into (4) leads to 
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Fig. 8. Filled and empty circles are hemispherical cloud fractions 
A, estimated by (14) and vertically projected cloud fractions + (1), 
respectively, plotted as functions of corresponding A, as determined 
directly from the Monte Carlo simulations for all 45 scenes. 


which converges for a 3 > — 2, though judging from Fig. 
7a, it appears unlikely that a 3 will be < —0.5. Regard- 
less of the fact that (13) yields A ( (/x) > 1, the small 
values of a 3 limit this usually to very small jx. This has 
little impact on estimates of A c , as can be seen in Fig. 
8, which compares A r determined by (14) with that ob- 
tained directly from the MC simulations. The agreement 
is almost perfect as all points lie virtually atop the 1:1 
line (despite a 3 being based on values for fi > 0.2). 
When a 3 is parameterized by the regression line (see 
Fig. 7a) 

a 3 « 0.17 lnA t .(l) (15) 

and used in (14), estimates of A c are almost as good as 
those shown in Fig. 8. Incidentally, use of (15) in (13) 
implies that, regardless of A c ( 1), A r ( /x) > 1 only for fx 
< £-i/o 17 ^ 0 Q028. Figure 8 also shows that the largest 
difference between A c and corresponding A,.(l) is —0.1 
(or 20%), which is for the same scene responsible for 
the anomalous values in Figs. 4, 5, and 7 [see scene 
Bll in Barker et al/s (1996) Fig. 1]. For the most part, 
however, relative differences between A c and A c ( 1) are 
typically only about 5%. Since errors in cloud fraction 
affect both clear and cloudy components of (5a), it may 
be necessary, at times, to consider distinguishing be- 
tween A ( . and A r (l) in GCMs. As shown later, something 
as simple as (14) and (15) will likely suffice (see ap- 
pendix A for an alternate approach). 

The final part of this section compares the magnitudes 
of all-sky transmittance biases due to (i) neglect of cloud 
sides [i.e., use of A t (l) as opposed to A ( ], and (ii) neglect 


of horizontal variable r (i.e., use of F|ft as opposed to 
TJE). This is achieved by examining the transmittances 
biases 

cloud side bias: 

A7V = {1 - A r ( 1) + A, (1)7^} - T 
variable r bias: 

A7; = {1 - A c + A ( T c Pf d (r)} - F, (16) 

where F is all-sky transmittance from the MC simula- 
tions and the terms in braces are approximate formulas. 
Hence, A T A is informed of true Monte Carlo cloud 
transmittances but has no information about cloud sides. 
On the other hand, AF T has Monte Carlo information 
about cloud sides but no information about horizontal 
variability of r. Figure 9 shows AT a and AF T for all 45 
scenes plotted against A r (l), r, and v. Clearly, the dom- 
inant bias is that due to variable r as the magnitude of 
A T t is typically 2-5 times larger than A T A for A ( .(l) ^ 

0. 9. The largest values of A T t are between -0.1 and 
—0.15 and occur for scenes with A r (l) near 0.75 (these 
are roughly 30%-50% relative biases, given that F for 
these cases are —0.4). The preference for A F T to be 
maximal for v ^ 1 is the result of a balance between 
having sufficiently many clouds to impact strongly the 
all-sky signal, but not too much cloud for, as A r (l) — » 

1, horizontal variability tends to weaken (Barker et al. 
1996). Likewise, despite high variability (small v) when 
A r (l) is very small (Barker et al. 1996), clouds con- 
tribute weakly to all-sky transmittance, thus reducing 
AF T . Relative biases for all-sky emissivity due to neglect 
of variable r are between +10% and +30% for the 
majority of scenes. 

Conversely, Fig. 9 shows that A T A tends to be greatest 
for A r (l) near 0.3, which often have v < 1. When vand 
A r (l) are small, and r is even just moderately large, 
there are sufficiently many deep clouds to initiate a large 
zenith angle dependence on cloud fraction (see Fig. 7a). 
Since this bears directly on the weighting of clear-sky 
transmittance, the cloud side bias is understandably larg- 
est for small cloud fractions. Of the 45 scenes, only 
scene Bll [see Barker et al.’s (1996) Fig. 1] has !AF a I 
> IAF T I (again the scene responsible for the anomalous 
points in previous plots). Having been viewed at large 
0 o , inferred values of r 0 , 83 can be anomalously large on 
the sunlit side of clouds (Barker and Liu 1995), and this 
would lead to excessive values of both h and a 3 . Also, 
note that while A7^ decreases slowly with increasing 
A r (l), vanishing for overcast, AF T for near overcast con- 
ditions ranges from 0 to about -0.06. 

Since AF^ and A T t are of opposite sign, neglect of 
both cloud sides and variable t (as in conventional PPH 
models) will yield biases between those plotted in Fig. 
9 but with a strong tendency to be negative (i.e., too 
little transmittance). Hence, it can be expected that 
GCMs under- and overestimate transmittances and em- 
issivities, respectively, for MBL cloud fields. The main 
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conclusion of this section, therefore, is that a simple 
parameterization for A ( , like (14) and (15) for example, 
should suffice, while more attention should be paid to 
the impact of horizontally variable r. This is addressed 
in the next section, which presents a simple parameter- 
ization for T c]d . 


5. A parameterized model for T cld 

Having established that (5a) and p(r II), hereinafter 
referred to as simply p(r), are likely to be adequate 
approximations for MBL clouds, this section presents a 
parameterization for T dd , as defined in (5b), that may 
be useful for climate modeling studies. It is essentially 
an independent pixel approximation IPA based on the 
assumption that frequency distributions of r often follow 
gamma distributions. 

Barker et al. (1996) demonstrated that for MBL 
clouds, Landsat-inferred p(f) are often represented very 
well by the normalized gamma distribution function, 
which can be written as 

r v x e " r ' T ; {r > 0; v > 0}, 

(17) 

where T(^) is the gamma function. Chambers et al. 
(1997) have demonstrated that Landsat-infened /?(r), as 
used by Barker et al. (1996), are reliable. Moreover, 
p(r ) produced by cloud resolving models for MBL 
clouds are also described well by p v (f) (S. Krueger and 
B. Stevens, 1996, personal communication). Addition- 
ally, Barker (1996) derived an IPA for computing solar 
radiative fluxes for horizontally inhomogeneous MBL 
clouds based on p v {f ). The parameter v defines the form 
of p, (r), but lacks a unique definition. For example. 
Barker et al. (1996) used the method of moments, as 
above, to estimate v as 

^ = {flay. (18a) 

Alternatively, the maximum likelihood estimate requires 
solving 


’ r>j 


where 


ilf(v) + In 



— lnr = 0, 


(18b) 


= T" InlV). (18c) 

dv 

While the value of ^ depends on the method of solution, 
differences are generally less than 20%. 2 


Fig. 9. Transmittance biases due to omission of (i) cloud side 
effects, A T A : and (ii) horizontally variable r. A T r plotted as functions 

of corresponding (a) A,(l), (b) r, and (c) v for all 45 scenes. See text 2 Exponentiating and rearranging (18b) leads to e lnT /r = e* ivt tv, 
and (16) for details. which also equals the r reduction factor in Cahalan et al.’s (1994a) 

effective thickness approximation. 
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Substituting (17) into (5b), it can be shown (see ap- 
pendix B) that 


7 , !' ld = jl - (1 - X) 

= 1 _ e cld» 

where 


v - (v + 1)t X 


;t' ,+ 


„=o v + 2 + nJJ 
(19) 


x 


V 

, — * 

V + T 


and the superscript Y indicates that (17) has been used. 
Note that the leading term on the right of (19) is trans- 
mittance for normal incident radiance and approaches 
e 1 as v — > oo. Also, as t — » 7*Jy ~ f~ v t whereas 
T? f d ~ e~Vf . For typical values of r and v, the series 
in (19) converges to 10“ 4 in less than 10 terms and often 
in less than five terms. This confines errors in T[ ld to 
the fifth decimal place, thus making accurate determi- 
nation of (19) efficient. As listed in appendix C, com- 
putational requirements of 7[ Jd are often 1-5 times those 
of conventional methods of computing r^p d . The only 
time (19) is overly cumbersome is when ris small and 
v is large: these conditions, however, appear to be rare 
(Barker et al. 1996). Moreover, when v is large it is 
adequate to revert to 7*£f d . When v is an integer, the 
infinite sum in (19) can be replaced by the finite sum 
(see appendix B) 


1 


x n+l 

V + 2 + n 


1 


jt" +l 


ln(l — jc) + 



( 20 ) 


which reduces greatly the time required to compute 
7^ ld (see appendix C). This identity could be used glob- 
ally if one is willing to round off vto the nearest integer. 
This is not necessarily as harsh as it sounds given the 
magnitude of other uncertainties in GCM cloud prop- 
erties. Rounding v < 1 to v = 1 would be undesirable 
given that dF[ ld Idv changes rapidly for v < 1 . Moreover, 
Table Cl shows that for v < 1, (19) is efficient. There- 
fore, the rounding need only be done for v > 1, where 
dT l clli ldv is relatively small. On the other hand, in an 
operational setting it may be desirable to use a 3D look- 
up table for (19). 

Figure 10 shows that 7*J?f d underestimates ( = 2 
djx) by about an order of magnitude for 7*J55 ^ 
0.2 and 20%-100% otherwise. In fact, the overall values 
of mean bias error (MBE) and root-mean-square error 
(rmse) are —0.100 and 0.131, respectively. Conversely, 
overall values of MBE and rmse for F[ ld relative to 
7*155 are +0.023 and 0.057, respectively. Thus, use of 
F[, d reduces the bias error by a factor of ~4 and the 
random error by a factor of ~2. Note that the positive 
bias for F[ ld increases as 7^5 increases [i.e., as v and r 
tend to decrease; see Fig. 9 and Barker et al. (1996)]. 
This is due, in part, to the fact that 7*[ ld is based on r 
e (0, oo) ? while 7^5 is based on Landsat-inferred r, which 




Fig. 10. Gamma-weighted transmittances T l M computed by (19) 
using r and v, and PPH transmittances TPf d computed using r for 45 
Landsat scenes plotted as functions of Monte Carlo transmittances 
The only difference between the plots is the scale of their axes. 
Solid lines represent perfect model performance, while broken lines 
are regression fits to the points. In (a), regression lines are of the 
form a + and are much determined by large transmittance cases. 
In (b), regression lines are of the form a (T^) h and are determined 
by small transmittance cases. 


has a minimum optical depth r mjn > 0. Thus, when v is 
small (<1), a significant contribution to F[ ld can come 
from 0 < r < r min . Wielicki and Parker (1992) estimate 
that the fractional amount of optically thin MBL cloud 
(r ^ 0.1) undetected by Landsat visible reflectance 
thresholds is typically less than 0.05. Since transmit- 
tances for r < r min are almost 1.0, the parameterized 
model interprets this as almost cloudless sky and so 
overestimates 7*$. Therefore, to call this strictly an over- 
estimation on the part of 7*[ ld is not entirely true. 

Figure 1 1 shows z£ ld and as defined in (19) and 
(6), plotted as functions of A c . Regardless of A t ., is 
almost always greater than 0.8. The stippled ellipse con- 
tains most of the Luo et al. (1994) estimates of daytime 
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Fig. 1 1. Mean cloud emittances predicted by the gamma-weighted 
PPH e' M and PPH e% models [see (19) and (6), respectively] as a 
function of hemispherical cloud fraction A t for 45 Landsat scenes. 
The stippled ellipse encompasses most of Luoet al.’s ( 1994) estimates 
of 1 1-/Ltm cloud emissivities for (250 km) 2 regions of ~l-km reso- 
lution AVHRR data viewed off the west coast of South America. 


11 -jLtm cloud emissivities for marine stratocumulus off 
the west coast of South America (see their Fig. 8). While 
their results are for (250 km) 2 images of ~ 1 km reso- 
lution Advanced Very High Resolution Radiometer 
data, they agree quite well with £r[ |d . Values of z£ ld for 
A c =£ 0.4 may, in reality, be shifted slightly down and 
right (i.e., closer to Luo et al.’s ellipse, which is not 
subject to the same shift as they used 1 1-pm radiances). 
This is because scenes with small cloud fractions (and 
small v) tend to have many small values of r and very 
thin clouds were certainly missed by Landsat (horizontal 
shift of probably <0.05) and as a result, r and v (and 
thus mean emittance) were overestimated slightly (ver- 
tical shift down). 

Finally, reconsider the prescription of cloud thickness 
h that was used to create 3D cloud fields [Eq. (8)]. 
Substituting (8) into (17) yields standard deviations of 
h, as shown in Fig. 12. Also shown are regions that 
contain the majority of scenes in three classes: (A) over- 
cast stratocumulus, (B) broken stratocumulus, and (C) 
scattered cumulus (see Table 3 in Barker et al. 1996). 
For the most part, standard deviations of h are between 
50 m and 125 m, which is in agreement with some 
observations (Loeb et al. 1997, manuscript submitted 
to J. Atmos. Sci.) and also fits well with a theoretical 
model (Considine et al. 1997). 


standard deviation of cloud geometric thickness (m) 



mean optical depth (1 1 .5 pm) 

Fig. 12. Standard deviation of cloud geometric thickness h (in 
meters) as a function of f and v. These were determined by consid- 
ering h to be related to r as in (8) and t to be defined by the gamma 
distribution in (17). Stippled regions demarcate where most of the 
scenes lie that were used by Barker et al. (1996): (A) overcast; (B) 
broken stratocumulus; (C) scattered cumulus. 

6. Summary and conclusions 

This paper presented a simple conceptual model for 
flux transmittance of longwave (LW) radiation 
through an inhomogeneous, marine boundary layer 
(MBL) cloud field. Two general aspects of cloud ge- 
ometry were addressed: horizontal variability of op- 
tical depth rand cloud sides. Using a 3D Monte Carlo 
photon transport algorithm and fields of r inferred 
from 45 Landsat images, it was demonstrated that 
when cloud fraction is ^0.9, neglect of horizontal 
variable r leads to all-sky transmittance biases that 
are roughly 2-5 times larger than, and opposite in 
sign to, biases stemming from neglect of cloud sides. 3 
As such, priority (in both research effort and com- 
putation) should be given to parameterizations that ac- 
count for the impact of variable r. 

Regarding horizontal variability of r, an approximate 
method for computing LW flux transmittances was fur- 
nished. It is essentially a stochastic radiative transfer 
model whose validity rests on the assumption that fre- 
quency distributions of optical depth p(r ), for (60 km) 2 
regions, are often approximated well by gamma distri- 


3 This disparity in biases may be much ameliorated, or even re- 
versed, for land cumulus, which exhibit both sharper edges (Wielicki 
and Parker 1992) and greater vertical extent than MBL clouds. This 
is the subject of a later study. 





2796 


JOURNAL OF THE ATMOSPHERIC SCIENCES 


Volume 54 


bution functions (Barker et al. 1996). This method is 
computationally efficient and suitable for GCMs. It was 
also demonstrated that the standard plane-parallel, ho- 
mogeneous (PPH) model often underestimates cloud 
transmittances by about an order of magnitude for thick 
clouds and by 20%-100% for thinner clouds. Converse- 
ly, when the mean and standard deviation of r are used 
to define a gamma distribution, the gamma-weighted 
PPH removes typically more than 80% of the homo- 
geneous bias. While this model neglected scattering, 
inclusion of it should not alter results much (a similar 
model with scattering is under development). 

It was shown that, in principle, cloud fraction and 
radiation cannot be decoupled. But, for computation of 
fluxes for MBL clouds, the conventional technique of 
weighting clear- and cloudy-sky transmittances by suit- 
able clear- and cloudy-sky fractions is acceptable (cf. 
Stephens 1988). But what is a suitable cloud fraction? 
Real clouds have depth and therefore, vertically pro- 
jected cloud fractions A c differ from cloud fractions A c 
presented to an isotropic beam of radiation. It is not 
obvious what cloud fractions GCM modelers think their 
radiation models are, and should, be using. Here A c is 
the more relevant quantity for computing fluxes, but it 
can be expected to be a complex function of cloud aspect 
ratios and spatial arrangement of clouds. Despite this, 
a simple parameterization of A ( as a function of A c was 
offered based on the 45 MBL cloud fields used here. 

Thus, the main recommendation stemming from this 
study is: the LW radiative effects of horizontal variable 
rfor MBL clouds should be included in GCM radiation 
routines in conjunction with due consideration of the 
effects of enhanced cloud fraction arising from hemi- 
spherical integration of cloud side view-factors. Since 
the combined effect of horizontal variability of r and 
cloud sides is to reduce cloud emittance relative to PPH 
conditions, the immediate impact of using the param- 
eterizations presented here in a GCM would be a slight 
warming of MBL clouds (which would not be undesir- 
able for many GCMs). As a final note, while the mag- 
nitude of LW biases presented here are comparable to 
their solar counterparts (Cahalan et al. 1994a; Barker et 
al. 1996), LW biases may dominate at times as they act 
continuously. 


A,(/x)- 1 - [1 -A t { \)}e 




= 1 — [1 — A t (1 )] exp 


Vi - /x 2 


-a. 




. (Al) 


where a 4 is a coefficient and 0 is zenith angle. This 
corresponds to identical cylinders distributed on a plane 
according to Poisson’s law (Avaste et al. 1974) and has 
been used to describe A ( (/x) by Ellingson (1982) and 
Barker et al. (1993), while Otterman (1984) used it in 
his vegetation albedo model to describe direct-beam in- 
terception. The reason this might be considered more 
attractive than (13) is because it confines A c {p) e [0, 
1]. Substituting (Al) into (4) yields 


A r ~ A,( 1) + [1 - A,( 1)] 


X a 4 [Ci(a 4 ) sin (a 4 ) - si (a 4 ) cos(a 4 )] 

= A c (l) + [1 -A c (\)]f(a 4 \ (A2) 


where Ci is the cosine integral and si = Si -ir!2 in 
which Si is the sine integral (Abramowitz and Stegun 
1964). For the 45 Landsat scenes, a 4 was fitted with 


a 4 ~ 0.03 + 0.07A ( .(1), (A3) 

and f(a 4 ) was parameterized by 


f(a 4 ) ~ 


-29.194a 4 

-20.248 - a 4 (21J29 + a 4 )’ 


a 4 € [0, 1], 
(A4) 


which has a maximum error of 0.0025. Hence, (A2) 
through (A4) is almost as efficient as (14) and (15). 

The reason why (13) through (15) was used in this 
study rather than this technique was simply because 
it performed slightly better. This is not to say that the 
method presented in this appendix performed poorly; 
it just had a minor tendency to overestimate Monte 
Carlo values of A t . Nonetheless, it was presented here 
anyway as future studies (such as with land cumulus, 
perhaps) might find this approach more appropriate. 


APPENDIX B 
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APPENDIX A 

An Alternate Approach to Estimate A c 

A somewhat more attractive technique for describing 
A t (/x), which fits the MC results about as well as (13), 
is 


Derivation of Eqs. (19) and (20) 

Substituting (17) into (5b) and carrying out the in- 
tegration with respect to jjl first yields 


2 (v 

TcId " nv)\r 


E,(T)T" 'e ” /T dr, (Bl) 


Substituting 


e A t ) = T - re 7 + r 2 £,(r)} (B2) 


into (Bl) gives 
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71 


1 / vY 17 T V 

r(p)\T/ 4 rj 

I 


i» — — r(v + i) 

V 4 T 


4- I " /r drj 

From Gradshteyn and Ryzhik (1980), 
£’ i (t)t , ' m £ ,,T/T dr 

T(v + 2) 


r 


p 4 r 


v 4 2 


x ->£, I 1 , ^ + 2, j' + 3; 


^ 4 t 


(B3) 


(B4) 


where 2 Ffa, b, c; x) is the hypergeometric function. 
Since 1 4 (*z 4 2) — (iz 4 3) = 0 in (B4), 2 F i in (B4) 
converges for all r > 0 and v > 0 (Arfken 1970). Also 
for 2 F y in (B4) 


I 1, v 4 2, v 4 3; 


v 4 r 


* O' + 2) 2 


„=(, v424w\y4T 


(B5) 


Substituting (B5) into (B4), and (B4) into (B3), and 
rearranging yields (19) in the text. 

When v is an integer, the infinite sum in (19) can be 
replaced by the finite sum in (20) for which a proof is 
given here. Let the sum in (19) be represented by 




-2 


" 0 n + (v + 2) ", m + (v + 1) 

Multiplying both sides of (B6) by x yfI gives 

^ x k 


x- k 'S = E 


m + {v 


— = E , 

1) <-4 + 2 k 


(B6) 


(B7) 


which can also be written as 


x" 


X 

■s-2 

*=1 




(B8) 


Since 


f 7 = — ln(l - x), 

t i k 

(B8) can be rewritten as 


(B9) 


Table Cl. CPU time required by an HP 725/75 workstation to 
compute T% via both (Cl) and (C2) normalized to the time required 
to compute (C2). Times for T l M are also normalized to those required 
for (C2). For T' M , values not in parentheses are for ( 19), while those 
in parentheses are for (20). For both Tjft [using Eq. (Cl)] and T[ u 
[using Eq. (19)], summations were terminated w'hen the magnitude 
of the terms became less than 0.0001. 


T 

[using 

(Cl)] 

T% 

[using 

(C2)l 


T[ u 


v = 0A 

v = 1.0 

v = 5.0 

0.46 

1.30 

1.0 

1.18 

3.76 (1.29) 

11.09(1.75) 

2.3 

0.66 

1.0 

1.18 

1.87 (1.29) 

4.06(1.75) 

6.9 

0.66 

1.0 

1.18 

1.52 (1.29) 

2.47 (1.75) 


APPENDIX C 

Computational Considerations for rj ld and 

This appendix documents some computational con- 
siderations of the homogeneous solution rj?f d which util- 
izes £ 3 (t), and the gamma- weighted solution T ] M as de- 
fined in (19). Here £ 3 (t) is evaluated using 

Ey(r) 


0.250621 4 2.334733t 4 t 2 


(1.681534 4 3.330657r 4 r 2 )re 7 


U- - y - lnr - X 


2 \2 


(-T)" 


i> (m - 2 )m\ ’ 


for r ^ 1 


for r < 1, 


(Cl) 


as given by (Abramowitz and Stegun 1964), and also 
by the parameterization 


E y (r) = exp 



for t € [0, 100], (C2) 


which has maximum errors of 0.2% and 1.37% for r 
< 55 and t < 100, respectively, and 


a , = 2277326.0 
a 2 = 18451521.0 
a, = 36528961.0 
a 4 = 15688104.0 
a, = 1115332.6 


b x = -3285487.0 
b 2 = -17343648.0 
by = -11795192.0 
b 4 = -1086098.3 
b 5 = -253.0 

K = 1.0. 


Table Cl lists some CPU requirements for compu- 
tation of and rj’ u . 


S 


1 


x 


v ± 1 


ln(l - x) 4 



(B10) 


which is equivalent to (20) and completes the proof. 
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